
library(tidyverse)
library(haven)
library(fixest)
library(splines)
library(ggsci)
library(ggthemes)
library(srvyr)
library(survey)
library(EnvStats)
library(arsenal)


# Data -----
set.seed(24)

india_pop_mortality <-
  read_xlsx(
    "~/BurkeLab Dropbox/Carlos Gould/In Praise/Data/WPP2022_GEN_F01_DEMOGRAPHIC_INDICATORS_COMPACT_REV1.xlsx",
    sheet = 2
  ) %>%
  filter(`Region, subregion, country or area *` == "India") %>%
  dplyr::select(
    `Region, subregion, country or area *`,
    Year,
    `Total Population, as of 1 January (thousands)`,
    `Total Deaths (thousands)`,
    `Crude Death Rate (deaths per 1,000 population)`
  ) %>%
  dplyr::rename(# year=Year,
    population = `Total Population, as of 1 January (thousands)`,
    deaths = `Total Deaths (thousands)`,
    death_rate = `Crude Death Rate (deaths per 1,000 population)`) %>%
  dplyr::select(Year, deaths, population, death_rate) %>%
  mutate_all(as.numeric) %>%
  mutate(deaths = deaths * 1000,
         population = population * 1000)

ires <- read_dta("~/BurkeLab Dropbox/Carlos Gould/Stanford/PMUY-AmbientAir/Data/IRES/dataverse_files (1)/CEEW - IRES Data.dta")


ires_predictions <- 
  ires %>% 
  mutate(
    
    yearly_cost_2019_actual_subsidized = as.numeric(q505_1_lpg_large_n)*(as.numeric(q508_lpg_last_refill_pay)-200),
    yearly_cost_2019 = as.numeric(q505_1_lpg_large_n)*750,
    yearly_cost_2019_subsidized = as.numeric(q505_1_lpg_large_n)*550,
    
    yearly_cost_2023 = as.numeric(q505_1_lpg_large_n)*1100,
    yearly_cost_2023_subsidized = as.numeric(q505_1_lpg_large_n)*900,
  ) %>% 
  mutate(
    percent_cost_2019 =  yearly_cost_2019/ (q234_month_exp*12),
    percent_cost_2019_subsidized = yearly_cost_2019_subsidized/ (q234_month_exp*12),
    
    percent_cost_2023 = yearly_cost_2023/ (q234_month_exp*12*(1.05224^4)),
    percent_cost_2023_subsidized = yearly_cost_2023_subsidized/ (q234_month_exp*12*(1.05224^4)),
    
    refills_2023_4_cost = 1100*4,
    refills_2023_9_cost = 1100*9
  ) %>% 
  mutate(
    refill_cat = ifelse(as.numeric(q505_1_lpg_large_n)<=3, "<4",
                        ifelse(as.numeric(q505_1_lpg_large_n)>3 & 
                                 as.numeric(q505_1_lpg_large_n<=9), "4-9",
                               ifelse(as.numeric(q505_1_lpg_large_n)>9, ">9", NA)))
  ) %>% 
  mutate(
    q504_lpg_pmuy_yn = ifelse(q504_lpg_pmuy_yn==99, NA, q504_lpg_pmuy_yn),
    
    yearly_cost_2019_actual_subsidized_2023usd = yearly_cost_2019_actual_subsidized*(1.05224^4),
    
    percent_cost_2019 = 750 / q234_month_exp,
    percent_cost_2019_subsidized= 550 / q234_month_exp,
    
    refill_cost_2023 = 1100,
    refill_cost_2023_subsidy = 1100-200,
    
    yearly_exp = 12*q234_month_exp,
    yearly_exp_2023 = 12*q234_month_exp*(1.05224^4),
    
    percent_cost_2023_4 = refills_2023_4_cost / (12*q234_month_exp*(1.05224^4)),
    percent_cost_2023_9 = refills_2023_9_cost / (12*q234_month_exp*(1.05224^4)),
    
    percent_cost_2023 = 1100 / (q234_month_exp*(1.05224^4)),
    percent_cost_2023_subsidized = 900 / (q234_month_exp*(1.05224^4)),
  ) %>% 
  mutate(
    refill_cat = factor(refill_cat, levels=c("<4", "4-9", ">9"))
  ) %>% 
  mutate(
    yearly_cost_2019_actual_subsidized_2023usd =
      ifelse(is.na(yearly_cost_2019_actual_subsidized_2023usd), (1.05224^4)*550, yearly_cost_2019_actual_subsidized_2023usd)
  )

ires_predictions_pmuy <-
  ires_predictions %>% 
  filter(q504_lpg_pmuy_yn==1) 

# weighted.mean(ires_predictions_pmuy$q213_no_members, w=ires_predictions_pmuy$sw_state)

# DEVELOP SCENARIOS --------

percent_price_change_1100 = (1100 - 550) / 550
percent_price_change_900 = (900 - 550) / 550
percent_price_change_700 = (700 - 550) / 550

# ESTIMATING LPG KG PER YEAR FROM PRICE ELASTICITIES --------

elasticity1 = rnorm(1, .33, .025)
elasticity2 = rnormTrunc(1, .2, .025, 0.1, elasticity1)
elasticity3 = rnormTrunc(1, .10, .0025, 0.01, elasticity2)

ires_predictions_pmuy1 <-
  ires_predictions_pmuy %>%
  mutate(kg_lpg_2019 = as.numeric(q505_1_lpg_large_n) * 14.2) %>%
  mutate(kg_lpg_2023_550 = as.numeric(q505_1_lpg_large_n) * 14.2) %>%
  mutate(
    kg_lpg_2023_1100 =
      ifelse(kg_lpg_2019 < 85.2,
        kg_lpg_2019 -
          (kg_lpg_2019 * (percent_price_change_1100 * elasticity1)),
        ifelse(
          kg_lpg_2019 >= 85.2 & kg_lpg_2019 < 127.8,
          kg_lpg_2019 -
            (kg_lpg_2019 * (percent_price_change_1100 * elasticity2)),
          ifelse(
            kg_lpg_2019 >= 127.8 & kg_lpg_2019 < 170.4,
            kg_lpg_2019 -
              (kg_lpg_2019 * (percent_price_change_1100 * elasticity3)),
            ifelse(kg_lpg_2019 >= 170.4, kg_lpg_2019, NA)
          )
        )
      ),
    kg_lpg_2023_900 =
      ifelse(kg_lpg_2019 < 85.2,
             kg_lpg_2019 -
               (kg_lpg_2019 * (percent_price_change_900 * elasticity1)),
             ifelse(
               kg_lpg_2019 >= 85.2 & kg_lpg_2019 < 127.8,
               kg_lpg_2019 -
                 (kg_lpg_2019 * (percent_price_change_900 * elasticity2)),
               ifelse(
                 kg_lpg_2019 >= 127.8 & kg_lpg_2019 < 170.4,
                 kg_lpg_2019 -
                   (kg_lpg_2019 * (percent_price_change_900 * elasticity3)),
                 ifelse(kg_lpg_2019 >= 170.4, kg_lpg_2019, NA)
               )
             )
      ),
    kg_lpg_2023_700 =
      ifelse(kg_lpg_2019 < 85.2,
             kg_lpg_2019 -
               (kg_lpg_2019 * (percent_price_change_700 * elasticity1)),
             ifelse(
               kg_lpg_2019 >= 85.2 & kg_lpg_2019 < 127.8,
               kg_lpg_2019 -
                 (kg_lpg_2019 * (percent_price_change_700 * elasticity2)),
               ifelse(
                 kg_lpg_2019 >= 127.8 & kg_lpg_2019 < 170.4,
                 kg_lpg_2019 -
                   (kg_lpg_2019 * (percent_price_change_700 * elasticity3)),
                 ifelse(kg_lpg_2019 >= 170.4, kg_lpg_2019, NA)
               )
             )
      )
  ) %>%
  mutate(
    kg_lpg_2023_1100 = ifelse(kg_lpg_2023_1100 < 0, 0, kg_lpg_2023_1100),
    kg_lpg_2023_900 = ifelse(kg_lpg_2023_900 < 0, 0, kg_lpg_2023_900),
    kg_lpg_2023_700 = ifelse(kg_lpg_2023_700 < 0, 0, kg_lpg_2023_700),
    kg_lpg_2023_550 = ifelse(kg_lpg_2023_550 < 0, 0, kg_lpg_2023_550)
  ) 

lpg_densities_fig <- 
  ggplot() + 
  geom_density(data=ires_predictions_pmuy1,
               aes(x=round(kg_lpg_2023_1100, digits=0)), adjust=1, bw=10) + 
  geom_density(data=ires_predictions_pmuy1,
               aes(x=round(kg_lpg_2023_900, digits=0)), adjust=1, bw=10, linetype="dashed") + 
  geom_density(data=ires_predictions_pmuy1,
               aes(x=round(kg_lpg_2023_700, digits=0)), linetype="dotted", adjust=1, bw=10) + 
  geom_density(data=ires_predictions_pmuy1,
               aes(x=round(kg_lpg_2023_550, digits=0)), linetype="twodash", adjust=1, bw=10) +
  scale_x_continuous(
    limits=c(0, 200),
    breaks=c(0, 14*2, 14*4, 14*6, 14*8, 14*10, 14*12, 14*14),
    labels=c("0\n0", "28\n2","56\n4","84\n6","112\n8","140\n10","168\n12", "196\n14")
  ) + 
  xlab("Yearly LPG consumption") + 
  ylab("density") + 
  # xlim(0, 200) + 
  # coord_cartesian(xlim=c(0, 200)) +
  theme_clean() +
  theme(plot.background = element_blank(),
        plot.title = element_blank())

  
specific_points <- seq(1, 200, 1)

pred_2023_1100 <- ggplot_build(lpg_densities_fig)$data[[1]]
pred_2023_900 <- ggplot_build(lpg_densities_fig)$data[[2]]
pred_2023_700 <- ggplot_build(lpg_densities_fig)$data[[3]]
pred_2023_550 <- ggplot_build(lpg_densities_fig)$data[[4]]

pred_2023_1100_density_values <- 
  approx(pred_2023_1100$x, pred_2023_1100$y, xout = specific_points)$y
pred_2023_900_density_values <- 
  approx(pred_2023_900$x, pred_2023_900$y, xout = specific_points)$y
pred_2023_700_density_values <- 
  approx(pred_2023_700$x, pred_2023_700$y, xout = specific_points)$y
pred_2023_550_density_values <-
  approx(pred_2023_550$x, pred_2023_550$y, xout = specific_points)$y

density_2019_2023 <- 
  bind_cols(
    specific_points,
    pred_2023_1100_density_values,
    pred_2023_900_density_values,
    pred_2023_700_density_values,
    pred_2023_550_density_values
  ) %>% 
  as_tibble() %>% 
  rename(
    kg_lpg = `...1`,
    pred_2023_1100_density_values = `...2`,
    pred_2023_900_density_values = `...3`,
    pred_2023_700_density_values = `...4`,
    pred_2023_550_density_values = `...5`
  ) %>% 
  mutate(
    sum_pred_2023_1100_density_values = sum(pred_2023_1100_density_values, na.rm=T),
    sum_pred_2023_900_density_values = sum(pred_2023_900_density_values, na.rm=T),
    sum_pred_2023_700_density_values = sum(pred_2023_700_density_values, na.rm=T),
    sum_pred_2023_550_density_values = sum(pred_2023_550_density_values, na.rm=T)
  ) %>% 
  mutate(
    ratio_sum_pred_2023_1100_density_values = 1/sum_pred_2023_1100_density_values,
    ratio_sum_pred_2023_900_density_values = 1/sum_pred_2023_900_density_values,
    ratio_sum_pred_2023_700_density_values = 1/sum_pred_2023_700_density_values,
    ratio_sum_pred_2023_550_density_values = 1/sum_pred_2023_550_density_values
  ) %>% 
  mutate(
    pred_2023_1100_density_values = pred_2023_1100_density_values * ratio_sum_pred_2023_1100_density_values,
    pred_2023_900_density_values = pred_2023_900_density_values * ratio_sum_pred_2023_900_density_values,
    pred_2023_700_density_values = pred_2023_700_density_values * ratio_sum_pred_2023_700_density_values,
    pred_2023_550_density_values = pred_2023_550_density_values * ratio_sum_pred_2023_550_density_values
  )


# Three PM2.5 exposure scenarios -------
min1 = rnorm(1, 25, 5)
max1 = rnorm(1, 100, 5)
decline1 = rnorm(1, 70, 2.5)
exp1 = rnorm(1, 4, .05)
exp2 = rnorm(1, 8, .1)

x=1:200
PM25_1=min1 + (max1-min1) / (1 + (x/decline1)^exp1)
PM25_2=min1 + (max1-min1) / (1 + (x/max1)^exp1)
PM25_3=max1 - ((max1 - min1) / 200)*x

kg_lpg_pm25 <- tibble(
  kg_lpg = 1:200,
  PM25_1 = PM25_1,
  PM25_2 = PM25_2,
  PM25_3 = PM25_3
)


# COMBINE KG LPG DENSITIES WITH PM25 SCENARIOS ------
india_kg_pm <- 
  density_2019_2023 %>% 
  left_join(kg_lpg_pm25)

# GENERATE FULL YEAR BY YEAR SCENARIOS --------
india_kg_pm_scenario <- 
  expand_grid(2023:2030,
              india_kg_pm) %>%
  rename(Year = `2023:2030`) %>%
  mutate(VSL = 640000,
         VSL_low = 330000,
         VSL_high = 1152800) %>%
  left_join(india_pop_mortality) %>%
  mutate(pmuy_fraction = (520.32 * 1000000) / population) %>%
  mutate(
    pred_1100_pop = population * pmuy_fraction * pred_2023_1100_density_values,
    pred_900_pop = population * pmuy_fraction * pred_2023_900_density_values,
    pred_700_pop = population * pmuy_fraction * pred_2023_700_density_values,
    pred_550_pop = population * pmuy_fraction * pred_2023_550_density_values
  )



# CALCULATE DAMAGES ------

india_kg_pm_scenario_1 <-
  india_kg_pm_scenario %>%
  mutate(PM25_1 = PM25_1 - 2.4,
         PM25_2 = PM25_2 - 2.4,
         PM25_3 = PM25_3 - 2.4) %>%
  mutate(rho1 = log(1 + (PM25_1 / 1.6)) / (1 + exp((15.5 - PM25_1) / 36.8)),
         rho2 = log(1 + (PM25_2 / 1.6)) / (1 + exp((15.5 - PM25_2) / 36.8)),
         rho3 = log(1 + (PM25_3 / 1.6)) / (1 + exp((15.5 - PM25_3) / 36.8))) %>%
  mutate(hazard1 = (0.1430) * rho1,
         hazard2 = (0.1430) * rho2,
         hazard3 = (0.1430) * rho3) %>%
  mutate(deathrate1 = (1 - (1 / exp(hazard1))) * death_rate,
         deathrate2 = (1 - (1 / exp(hazard2))) * death_rate,
         deathrate3 = (1 - (1 / exp(hazard3))) * death_rate) %>%
  mutate(
    pred_1100_deaths1 = (deathrate1 * pred_1100_pop) / 1000,
    pred_900_deaths1 = (deathrate1 * pred_900_pop) / 1000,
    pred_700_deaths1 = (deathrate1 * pred_700_pop) / 1000,
    pred_550_deaths1 = (deathrate1 * pred_550_pop) / 1000,
    
    pred_1100_deaths2 = (deathrate2 * pred_1100_pop) / 1000,
    pred_900_deaths2 = (deathrate2 * pred_900_pop) / 1000,
    pred_700_deaths2 = (deathrate2 * pred_700_pop) / 1000,
    pred_550_deaths2 = (deathrate2 * pred_550_pop) / 1000,
    
    pred_1100_deaths3 = (deathrate3 * pred_1100_pop) / 1000,
    pred_900_deaths3 = (deathrate3 * pred_900_pop) / 1000,
    pred_700_deaths3 = (deathrate3 * pred_700_pop) / 1000,
    pred_550_deaths3 = (deathrate3 * pred_550_pop) / 1000
    
  ) 

india_kg_pm_scenario_1_yr <-
  india_kg_pm_scenario_1 %>%
  mutate(
    pred_1100_kg = (pred_1100_pop * kg_lpg) / 5.5, # number of people at each kg lpg increment times kg lpg divided by number of people in each household
    pred_900_kg = (pred_900_pop * kg_lpg) / 5.5,
    pred_700_kg = (pred_700_pop * kg_lpg) / 5.5,
    pred_550_kg = (pred_550_pop * kg_lpg) / 5.5
  ) %>% 
  fgroup_by(Year) %>% 
  fsummarise(
    pred_1100_pop = fsum(pred_1100_pop, na.rm=T),
    pred_900_pop = fsum(pred_900_pop, na.rm=T),
    pred_700_pop = fsum(pred_700_pop, na.rm=T),
    pred_550_pop = fsum(pred_550_pop, na.rm=T),
    
    pred_1100_kg = fsum(pred_1100_kg, na.rm=T),
    pred_900_kg = fsum(pred_900_kg, na.rm=T),
    pred_700_kg = fsum(pred_700_kg, na.rm=T),
    pred_550_kg = fsum(pred_550_kg, na.rm=T),
    
    PM25_1_1100 = fmean(PM25_1, na.rm=T, w=pred_1100_pop),
    PM25_2_1100 = fmean(PM25_2, na.rm=T, w=pred_1100_pop),
    PM25_3_1100 = fmean(PM25_3, na.rm=T, w=pred_1100_pop),
    
    PM25_1_900 = fmean(PM25_1, na.rm=T, w=pred_900_pop),
    PM25_2_900 = fmean(PM25_2, na.rm=T, w=pred_900_pop),
    PM25_3_900 = fmean(PM25_3, na.rm=T, w=pred_900_pop),
    
    PM25_1_700 = fmean(PM25_1, na.rm=T, w=pred_700_pop),
    PM25_2_700 = fmean(PM25_2, na.rm=T, w=pred_700_pop),
    PM25_3_700 = fmean(PM25_3, na.rm=T, w=pred_700_pop),
    
    PM25_1_550 = fmean(PM25_1, na.rm=T, w=pred_550_pop),
    PM25_2_550 = fmean(PM25_2, na.rm=T, w=pred_550_pop),
    PM25_3_550 = fmean(PM25_3, na.rm=T, w=pred_550_pop),
    
    VSL = fmin(VSL, na.rm=T),
    VSL_low = fmin(VSL_low, na.rm=T),
    VSL_high = fmin(VSL_high, na.rm=T),
    
    pred_1100_deaths1 = fsum(pred_1100_deaths1, na.rm=T),
    pred_900_deaths1 = fsum(pred_900_deaths1, na.rm=T),
    pred_700_deaths1 = fsum(pred_700_deaths1, na.rm=T),
    pred_550_deaths1 = fsum(pred_550_deaths1, na.rm=T),
    
    pred_1100_deaths2 = fsum(pred_1100_deaths2, na.rm=T),
    pred_900_deaths2 = fsum(pred_900_deaths2, na.rm=T),
    pred_700_deaths2 = fsum(pred_700_deaths2, na.rm=T),
    pred_550_deaths2 = fsum(pred_550_deaths2, na.rm=T),
    
    pred_1100_deaths3 = fsum(pred_1100_deaths3, na.rm=T),
    pred_900_deaths3 = fsum(pred_900_deaths3, na.rm=T),
    pred_700_deaths3 = fsum(pred_700_deaths3, na.rm=T),
    pred_550_deaths3 = fsum(pred_550_deaths3, na.rm=T)

  ) %>%
  mutate(
    pred_900_deaths_difference1 = pred_1100_deaths1 - pred_900_deaths1,
    pred_700_deaths_difference1 = pred_1100_deaths1 - pred_700_deaths1,
    pred_550_deaths_difference1 = pred_1100_deaths1 - pred_550_deaths1,
    
    pred_900_deaths_difference2 = pred_1100_deaths2 - pred_900_deaths2,
    pred_700_deaths_difference2 = pred_1100_deaths2 - pred_700_deaths2,
    pred_550_deaths_difference2 = pred_1100_deaths2 - pred_550_deaths2,
    
    pred_900_deaths_difference3 = pred_1100_deaths3 - pred_900_deaths3,
    pred_700_deaths_difference3 = pred_1100_deaths3 - pred_700_deaths3,
    pred_550_deaths_difference3 = pred_1100_deaths3 - pred_550_deaths3
  ) %>%
  mutate(
    pred_900_deaths_difference_VSL1 = pred_900_deaths_difference1 * VSL,
    pred_900_deaths_difference_VSL1_low = pred_900_deaths_difference1 * VSL_low,
    pred_900_deaths_difference_VSL1_high = pred_900_deaths_difference1 * VSL_high,
    
    pred_900_deaths_difference_VSL2 = pred_900_deaths_difference2 * VSL,
    pred_900_deaths_difference_VSL2_low = pred_900_deaths_difference2 * VSL_low,
    pred_900_deaths_difference_VSL2_high = pred_900_deaths_difference2 * VSL_high,
    
    pred_900_deaths_difference_VSL3 = pred_900_deaths_difference3 * VSL,
    pred_900_deaths_difference_VSL3_low = pred_900_deaths_difference3 * VSL_low,
    pred_900_deaths_difference_VSL3_high = pred_900_deaths_difference3 * VSL_high,
    
    pred_700_deaths_difference_VSL1 = pred_700_deaths_difference1 * VSL,
    pred_700_deaths_difference_VSL1_low = pred_700_deaths_difference1 * VSL_low,
    pred_700_deaths_difference_VSL1_high = pred_700_deaths_difference1 * VSL_high,
    
    pred_700_deaths_difference_VSL2 = pred_700_deaths_difference2 * VSL,
    pred_700_deaths_difference_VSL2_low = pred_700_deaths_difference2 * VSL_low,
    pred_700_deaths_difference_VSL2_high = pred_700_deaths_difference2 * VSL_high,
    
    pred_700_deaths_difference_VSL3 = pred_700_deaths_difference3 * VSL,
    pred_700_deaths_difference_VSL3_low = pred_700_deaths_difference3 * VSL_low,
    pred_700_deaths_difference_VSL3_high = pred_700_deaths_difference3 * VSL_high,
    
    pred_550_deaths_difference_VSL1 = pred_550_deaths_difference1 * VSL,
    pred_550_deaths_difference_VSL1_low = pred_550_deaths_difference1 * VSL_low,
    pred_550_deaths_difference_VSL1_high = pred_550_deaths_difference1 * VSL_high,
    
    pred_550_deaths_difference_VSL2 = pred_550_deaths_difference2 * VSL,
    pred_550_deaths_difference_VSL2_low = pred_550_deaths_difference2 * VSL_low,
    pred_550_deaths_difference_VSL2_high = pred_550_deaths_difference2 * VSL_high,
    
    pred_550_deaths_difference_VSL3 = pred_550_deaths_difference3 * VSL,
    pred_550_deaths_difference_VSL3_low = pred_550_deaths_difference3 * VSL_low,
    pred_550_deaths_difference_VSL3_high = pred_550_deaths_difference3 * VSL_high
  ) %>%
  mutate(
    pred_900_deaths_difference_VSL1_discounted = pred_900_deaths_difference_VSL1 / (1 + 0.09) ^
      (Year - 2023),
    pred_900_deaths_difference_VSL1_low_discounted =  pred_900_deaths_difference_VSL1_low / (1 + 0.09) ^
      (Year - 2023),
    pred_900_deaths_difference_VSL1_high_discounted =  pred_900_deaths_difference_VSL1_high / (1 + 0.09) ^
      (Year - 2023),
    pred_900_deaths_difference_VSL2_discounted = pred_900_deaths_difference_VSL2 / (1 + 0.09) ^
      (Year - 2023),
    pred_900_deaths_difference_VSL2_low_discounted =  pred_900_deaths_difference_VSL2_low / (1 + 0.09) ^
      (Year - 2023),
    pred_900_deaths_difference_VSL2_high_discounted =  pred_900_deaths_difference_VSL2_high / (1 + 0.09) ^
      (Year - 2023),
    pred_900_deaths_difference_VSL3_discounted = pred_900_deaths_difference_VSL3 / (1 + 0.09) ^
      (Year - 2023),
    pred_900_deaths_difference_VSL3_low_discounted =  pred_900_deaths_difference_VSL3_low / (1 + 0.09) ^
      (Year - 2023),
    pred_900_deaths_difference_VSL3_high_discounted =  pred_900_deaths_difference_VSL3_high / (1 + 0.09) ^
      (Year - 2023),
    
    
    pred_700_deaths_difference_VSL1_discounted = pred_700_deaths_difference_VSL1 / (1 + 0.09) ^
      (Year - 2023),
    pred_700_deaths_difference_VSL1_low_discounted =  pred_700_deaths_difference_VSL1_low / (1 + 0.09) ^
      (Year - 2023),
    pred_700_deaths_difference_VSL1_high_discounted =  pred_700_deaths_difference_VSL1_high / (1 + 0.09) ^
      (Year - 2023),
    pred_700_deaths_difference_VSL2_discounted = pred_700_deaths_difference_VSL2 / (1 + 0.09) ^
      (Year - 2023),
    pred_700_deaths_difference_VSL2_low_discounted =  pred_700_deaths_difference_VSL2_low / (1 + 0.09) ^
      (Year - 2023),
    pred_700_deaths_difference_VSL2_high_discounted =  pred_700_deaths_difference_VSL2_high / (1 + 0.09) ^
      (Year - 2023),
    pred_700_deaths_difference_VSL3_discounted = pred_700_deaths_difference_VSL3 / (1 + 0.09) ^
      (Year - 2023),
    pred_700_deaths_difference_VSL3_low_discounted =  pred_700_deaths_difference_VSL3_low / (1 + 0.09) ^
      (Year - 2023),
    pred_700_deaths_difference_VSL3_high_discounted =  pred_700_deaths_difference_VSL3_high / (1 + 0.09) ^
      (Year - 2023),
    
    
    pred_550_deaths_difference_VSL1_discounted = pred_550_deaths_difference_VSL1 / (1 + 0.09) ^
      (Year - 2023),
    pred_550_deaths_difference_VSL1_low_discounted =  pred_550_deaths_difference_VSL1_low / (1 + 0.09) ^
      (Year - 2023),
    pred_550_deaths_difference_VSL1_high_discounted =  pred_550_deaths_difference_VSL1_high / (1 + 0.09) ^
      (Year - 2023),
    pred_550_deaths_difference_VSL2_discounted = pred_550_deaths_difference_VSL2 / (1 + 0.09) ^
      (Year - 2023),
    pred_550_deaths_difference_VSL2_low_discounted =  pred_550_deaths_difference_VSL2_low / (1 + 0.09) ^
      (Year - 2023),
    pred_550_deaths_difference_VSL2_high_discounted =  pred_550_deaths_difference_VSL2_high / (1 + 0.09) ^
      (Year - 2023),
    pred_550_deaths_difference_VSL3_discounted = pred_550_deaths_difference_VSL3 / (1 + 0.09) ^
      (Year - 2023),
    pred_550_deaths_difference_VSL3_low_discounted =  pred_550_deaths_difference_VSL3_low / (1 + 0.09) ^
      (Year - 2023),
    pred_550_deaths_difference_VSL3_high_discounted =  pred_550_deaths_difference_VSL3_high / (1 + 0.09) ^
      (Year - 2023)
    
  )


india_kg_pm_scenario_1_yr_long <- india_kg_pm_scenario_1_yr

# RUN BOOTS ------

for (i in 1:1000) {
  
  # ESTIMATING LPG KG PER YEAR FROM PRICE ELASTICITIES --------
  
  elasticity1 = rnorm(1, .33, .025)
  elasticity2 = rnormTrunc(1, .2, .025, 0.1, elasticity1)
  elasticity3 = rnormTrunc(1, .10, .0025, 0.01, elasticity2)
  
  ires_predictions_pmuy1 <-
    ires_predictions_pmuy %>%
    mutate(kg_lpg_2019 = as.numeric(q505_1_lpg_large_n) * 14.2) %>%
    mutate(kg_lpg_2023_550 = as.numeric(q505_1_lpg_large_n) * 14.2) %>%
    mutate(
      kg_lpg_2023_1100 =
        ifelse(kg_lpg_2019 < 85.2,
               kg_lpg_2019 -
                 (kg_lpg_2019 * (percent_price_change_1100 * elasticity1)),
               ifelse(
                 kg_lpg_2019 >= 85.2 & kg_lpg_2019 < 127.8,
                 kg_lpg_2019 -
                   (kg_lpg_2019 * (percent_price_change_1100 * elasticity2)),
                 ifelse(
                   kg_lpg_2019 >= 127.8 & kg_lpg_2019 < 170.4,
                   kg_lpg_2019 -
                     (kg_lpg_2019 * (percent_price_change_1100 * elasticity3)),
                   ifelse(kg_lpg_2019 >= 170.4, kg_lpg_2019, NA)
                 )
               )
        ),
      kg_lpg_2023_900 =
        ifelse(kg_lpg_2019 < 85.2,
               kg_lpg_2019 -
                 (kg_lpg_2019 * (percent_price_change_900 * elasticity1)),
               ifelse(
                 kg_lpg_2019 >= 85.2 & kg_lpg_2019 < 127.8,
                 kg_lpg_2019 -
                   (kg_lpg_2019 * (percent_price_change_900 * elasticity2)),
                 ifelse(
                   kg_lpg_2019 >= 127.8 & kg_lpg_2019 < 170.4,
                   kg_lpg_2019 -
                     (kg_lpg_2019 * (percent_price_change_900 * elasticity3)),
                   ifelse(kg_lpg_2019 >= 170.4, kg_lpg_2019, NA)
                 )
               )
        ),
      kg_lpg_2023_700 =
        ifelse(kg_lpg_2019 < 85.2,
               kg_lpg_2019 -
                 (kg_lpg_2019 * (percent_price_change_700 * elasticity1)),
               ifelse(
                 kg_lpg_2019 >= 85.2 & kg_lpg_2019 < 127.8,
                 kg_lpg_2019 -
                   (kg_lpg_2019 * (percent_price_change_700 * elasticity2)),
                 ifelse(
                   kg_lpg_2019 >= 127.8 & kg_lpg_2019 < 170.4,
                   kg_lpg_2019 -
                     (kg_lpg_2019 * (percent_price_change_700 * elasticity3)),
                   ifelse(kg_lpg_2019 >= 170.4, kg_lpg_2019, NA)
                 )
               )
        )
    ) %>%
    mutate(
      kg_lpg_2023_1100 = ifelse(kg_lpg_2023_1100 < 0, 0, kg_lpg_2023_1100),
      kg_lpg_2023_900 = ifelse(kg_lpg_2023_900 < 0, 0, kg_lpg_2023_900),
      kg_lpg_2023_700 = ifelse(kg_lpg_2023_700 < 0, 0, kg_lpg_2023_700),
      kg_lpg_2023_550 = ifelse(kg_lpg_2023_550 < 0, 0, kg_lpg_2023_550)
    ) 
  
  lpg_densities_fig <- 
    ggplot() + 
    geom_density(data=ires_predictions_pmuy1,
                 aes(x=round(kg_lpg_2023_1100, digits=0)), adjust=1, bw=10) + 
    geom_density(data=ires_predictions_pmuy1,
                 aes(x=round(kg_lpg_2023_900, digits=0)), adjust=1, bw=10, linetype="dashed") + 
    geom_density(data=ires_predictions_pmuy1,
                 aes(x=round(kg_lpg_2023_700, digits=0)), linetype="dotted", adjust=1, bw=10) + 
    geom_density(data=ires_predictions_pmuy1,
                 aes(x=round(kg_lpg_2023_550, digits=0)), linetype="twodash", adjust=1, bw=10) +
    scale_x_continuous(
      limits=c(0, 200),
      breaks=c(0, 14*2, 14*4, 14*6, 14*8, 14*10, 14*12, 14*14),
      labels=c("0\n0", "28\n2","56\n4","84\n6","112\n8","140\n10","168\n12", "196\n14")
    ) + 
    xlab("Yearly LPG consumption") + 
    ylab("density") + 
    # xlim(0, 200) + 
    # coord_cartesian(xlim=c(0, 200)) +
    theme_clean() +
    theme(plot.background = element_blank(),
          plot.title = element_blank())
  
  
  specific_points <- seq(1, 200, 1)
  
  pred_2023_1100 <- ggplot_build(lpg_densities_fig)$data[[1]]
  pred_2023_900 <- ggplot_build(lpg_densities_fig)$data[[2]]
  pred_2023_700 <- ggplot_build(lpg_densities_fig)$data[[3]]
  pred_2023_550 <- ggplot_build(lpg_densities_fig)$data[[4]]
  
  pred_2023_1100_density_values <- 
    approx(pred_2023_1100$x, pred_2023_1100$y, xout = specific_points)$y
  pred_2023_900_density_values <- 
    approx(pred_2023_900$x, pred_2023_900$y, xout = specific_points)$y
  pred_2023_700_density_values <- 
    approx(pred_2023_700$x, pred_2023_700$y, xout = specific_points)$y
  pred_2023_550_density_values <-
    approx(pred_2023_550$x, pred_2023_550$y, xout = specific_points)$y
  
  density_2019_2023 <- 
    bind_cols(
      specific_points,
      pred_2023_1100_density_values,
      pred_2023_900_density_values,
      pred_2023_700_density_values,
      pred_2023_550_density_values
    ) %>% 
    as_tibble() %>% 
    rename(
      kg_lpg = `...1`,
      pred_2023_1100_density_values = `...2`,
      pred_2023_900_density_values = `...3`,
      pred_2023_700_density_values = `...4`,
      pred_2023_550_density_values = `...5`
    )  %>% 
    mutate(
      sum_pred_2023_1100_density_values = sum(pred_2023_1100_density_values, na.rm=T),
      sum_pred_2023_900_density_values = sum(pred_2023_900_density_values, na.rm=T),
      sum_pred_2023_700_density_values = sum(pred_2023_700_density_values, na.rm=T),
      sum_pred_2023_550_density_values = sum(pred_2023_550_density_values, na.rm=T)
    ) %>% 
    mutate(
      ratio_sum_pred_2023_1100_density_values = 1/sum_pred_2023_1100_density_values,
      ratio_sum_pred_2023_900_density_values = 1/sum_pred_2023_900_density_values,
      ratio_sum_pred_2023_700_density_values = 1/sum_pred_2023_700_density_values,
      ratio_sum_pred_2023_550_density_values = 1/sum_pred_2023_550_density_values
    ) %>% 
    mutate(
      pred_2023_1100_density_values = pred_2023_1100_density_values * ratio_sum_pred_2023_1100_density_values,
      pred_2023_900_density_values = pred_2023_900_density_values * ratio_sum_pred_2023_900_density_values,
      pred_2023_700_density_values = pred_2023_700_density_values * ratio_sum_pred_2023_700_density_values,
      pred_2023_550_density_values = pred_2023_550_density_values * ratio_sum_pred_2023_550_density_values
    )
  
  # write_rds(density_2019_2023, "~/Downloads/india_subsidies_kg_density_2019_2023.rds")
  
  # Three PM2.5 exposure scenarios -------
  min1 = rnorm(1, 25, 5)
  max1 = rnorm(1, 100, 5)
  decline1 = rnorm(1, 70, 2.5)
  exp1 = rnorm(1, 4, .05)
  exp2 = rnorm(1, 8, .1)
  
  x=1:200
  PM25_1=min1 + (max1-min1) / (1 + (x/decline1)^exp1)
  PM25_2=min1 + (max1-min1) / (1 + (x/max1)^exp1)
  PM25_3=max1 - ((max1 - min1) / 200)*x
  
  kg_lpg_pm25 <- tibble(
    kg_lpg = 1:200,
    PM25_1 = PM25_1,
    PM25_2 = PM25_2,
    PM25_3 = PM25_3
  )
  

  
  # COMBINE KG LPG DENSITIES WITH PM25 SCENARIOS ------
  india_kg_pm <- 
    density_2019_2023 %>% 
    left_join(kg_lpg_pm25)
  
  # GENERATE FULL YEAR BY YEAR SCENARIOS --------
  india_kg_pm_scenario <- 
    expand_grid(2023:2030,
                india_kg_pm) %>%
    rename(Year = `2023:2030`) %>%
    mutate(VSL = 640000,
           VSL_low = 330000,
           VSL_high = 1152800) %>%
    left_join(india_pop_mortality) %>%
    mutate(pmuy_fraction = (520.32 * 1000000) / population) %>%
    mutate(
      pred_1100_pop = population * pmuy_fraction * pred_2023_1100_density_values,
      pred_900_pop = population * pmuy_fraction * pred_2023_900_density_values,
      pred_700_pop = population * pmuy_fraction * pred_2023_700_density_values,
      pred_550_pop = population * pmuy_fraction * pred_2023_550_density_values
    )
  
  
  
  # CALCULATE DAMAGES ------
  
  india_kg_pm_scenario_1 <-
    india_kg_pm_scenario %>%
    mutate(PM25_1 = PM25_1 - 2.4,
           PM25_2 = PM25_2 - 2.4,
           PM25_3 = PM25_3 - 2.4) %>%
    mutate(rho1 = log(1 + (PM25_1 / 1.6)) / (1 + exp((15.5 - PM25_1) / 36.8)),
           rho2 = log(1 + (PM25_2 / 1.6)) / (1 + exp((15.5 - PM25_2) / 36.8)),
           rho3 = log(1 + (PM25_3 / 1.6)) / (1 + exp((15.5 - PM25_3) / 36.8))) %>%
    mutate(hazard1 = (0.1430) * rho1,
           hazard2 = (0.1430) * rho2,
           hazard3 = (0.1430) * rho3) %>%
    mutate(deathrate1 = (1 - (1 / exp(hazard1))) * death_rate,
           deathrate2 = (1 - (1 / exp(hazard2))) * death_rate,
           deathrate3 = (1 - (1 / exp(hazard3))) * death_rate) %>%
    mutate(
      pred_1100_deaths1 = (deathrate1 * pred_1100_pop) / 1000,
      pred_900_deaths1 = (deathrate1 * pred_900_pop) / 1000,
      pred_700_deaths1 = (deathrate1 * pred_700_pop) / 1000,
      pred_550_deaths1 = (deathrate1 * pred_550_pop) / 1000,
      
      pred_1100_deaths2 = (deathrate2 * pred_1100_pop) / 1000,
      pred_900_deaths2 = (deathrate2 * pred_900_pop) / 1000,
      pred_700_deaths2 = (deathrate2 * pred_700_pop) / 1000,
      pred_550_deaths2 = (deathrate2 * pred_550_pop) / 1000,
      
      pred_1100_deaths3 = (deathrate3 * pred_1100_pop) / 1000,
      pred_900_deaths3 = (deathrate3 * pred_900_pop) / 1000,
      pred_700_deaths3 = (deathrate3 * pred_700_pop) / 1000,
      pred_550_deaths3 = (deathrate3 * pred_550_pop) / 1000
      
    ) 
  
  india_kg_pm_scenario_1_yr <-
    india_kg_pm_scenario_1 %>%
    mutate(
      pred_1100_kg = (pred_1100_pop * kg_lpg) / 5.5, # number of people at each kg lpg increment times kg lpg divided by number of people in each household
      pred_900_kg = (pred_900_pop * kg_lpg) / 5.5,
      pred_700_kg = (pred_700_pop * kg_lpg) / 5.5,
      pred_550_kg = (pred_550_pop * kg_lpg) / 5.5
      ) %>% 
    fgroup_by(Year) %>% 
    fsummarise(
      pred_1100_pop = fsum(pred_1100_pop, na.rm=T),
      pred_900_pop = fsum(pred_900_pop, na.rm=T),
      pred_700_pop = fsum(pred_700_pop, na.rm=T),
      pred_550_pop = fsum(pred_550_pop, na.rm=T),
      
      pred_1100_kg = fsum(pred_1100_kg, na.rm=T),
      pred_900_kg = fsum(pred_900_kg, na.rm=T),
      pred_700_kg = fsum(pred_700_kg, na.rm=T),
      pred_550_kg = fsum(pred_550_kg, na.rm=T),
      
      PM25_1_1100 = fmean(PM25_1, na.rm=T, w=pred_2023_1100_density_values),
      PM25_2_1100 = fmean(PM25_2, na.rm=T, w=pred_2023_1100_density_values),
      PM25_3_1100 = fmean(PM25_3, na.rm=T, w=pred_2023_1100_density_values),
      
      PM25_1_900 = fmean(PM25_1, na.rm=T, w=pred_2023_900_density_values),
      PM25_2_900 = fmean(PM25_2, na.rm=T, w=pred_2023_900_density_values),
      PM25_3_900 = fmean(PM25_3, na.rm=T, w=pred_2023_900_density_values),
      
      PM25_1_700 = fmean(PM25_1, na.rm=T, w=pred_2023_700_density_values),
      PM25_2_700 = fmean(PM25_2, na.rm=T, w=pred_2023_700_density_values),
      PM25_3_700 = fmean(PM25_3, na.rm=T, w=pred_2023_700_density_values),
      
      PM25_1_550 = fmean(PM25_1, na.rm=T, w=pred_2023_550_density_values),
      PM25_2_550 = fmean(PM25_2, na.rm=T, w=pred_2023_550_density_values),
      PM25_3_550 = fmean(PM25_3, na.rm=T, w=pred_2023_550_density_values),
      
      VSL = fmin(VSL, na.rm=T),
      VSL_low = fmin(VSL_low, na.rm=T),
      VSL_high = fmin(VSL_high, na.rm=T),
      
      pred_1100_deaths1 = fsum(pred_1100_deaths1, na.rm=T),
      pred_900_deaths1 = fsum(pred_900_deaths1, na.rm=T),
      pred_700_deaths1 = fsum(pred_700_deaths1, na.rm=T),
      pred_550_deaths1 = fsum(pred_550_deaths1, na.rm=T),
      
      pred_1100_deaths2 = fsum(pred_1100_deaths2, na.rm=T),
      pred_900_deaths2 = fsum(pred_900_deaths2, na.rm=T),
      pred_700_deaths2 = fsum(pred_700_deaths2, na.rm=T),
      pred_550_deaths2 = fsum(pred_550_deaths2, na.rm=T),
      
      pred_1100_deaths3 = fsum(pred_1100_deaths3, na.rm=T),
      pred_900_deaths3 = fsum(pred_900_deaths3, na.rm=T),
      pred_700_deaths3 = fsum(pred_700_deaths3, na.rm=T),
      pred_550_deaths3 = fsum(pred_550_deaths3, na.rm=T)
      
    ) %>%
    mutate(
      pred_900_deaths_difference1 = pred_1100_deaths1 - pred_900_deaths1,
      pred_700_deaths_difference1 = pred_1100_deaths1 - pred_700_deaths1,
      pred_550_deaths_difference1 = pred_1100_deaths1 - pred_550_deaths1,
      
      pred_900_deaths_difference2 = pred_1100_deaths2 - pred_900_deaths2,
      pred_700_deaths_difference2 = pred_1100_deaths2 - pred_700_deaths2,
      pred_550_deaths_difference2 = pred_1100_deaths2 - pred_550_deaths2,
      
      pred_900_deaths_difference3 = pred_1100_deaths3 - pred_900_deaths3,
      pred_700_deaths_difference3 = pred_1100_deaths3 - pred_700_deaths3,
      pred_550_deaths_difference3 = pred_1100_deaths3 - pred_550_deaths3
    ) %>%
    mutate(
      pred_900_deaths_difference_VSL1 = pred_900_deaths_difference1 * VSL,
      pred_900_deaths_difference_VSL1_low = pred_900_deaths_difference1 * VSL_low,
      pred_900_deaths_difference_VSL1_high = pred_900_deaths_difference1 * VSL_high,
      
      pred_900_deaths_difference_VSL2 = pred_900_deaths_difference2 * VSL,
      pred_900_deaths_difference_VSL2_low = pred_900_deaths_difference2 * VSL_low,
      pred_900_deaths_difference_VSL2_high = pred_900_deaths_difference2 * VSL_high,
      
      pred_900_deaths_difference_VSL3 = pred_900_deaths_difference3 * VSL,
      pred_900_deaths_difference_VSL3_low = pred_900_deaths_difference3 * VSL_low,
      pred_900_deaths_difference_VSL3_high = pred_900_deaths_difference3 * VSL_high,
      
      pred_700_deaths_difference_VSL1 = pred_700_deaths_difference1 * VSL,
      pred_700_deaths_difference_VSL1_low = pred_700_deaths_difference1 * VSL_low,
      pred_700_deaths_difference_VSL1_high = pred_700_deaths_difference1 * VSL_high,
      
      pred_700_deaths_difference_VSL2 = pred_700_deaths_difference2 * VSL,
      pred_700_deaths_difference_VSL2_low = pred_700_deaths_difference2 * VSL_low,
      pred_700_deaths_difference_VSL2_high = pred_700_deaths_difference2 * VSL_high,
      
      pred_700_deaths_difference_VSL3 = pred_700_deaths_difference3 * VSL,
      pred_700_deaths_difference_VSL3_low = pred_700_deaths_difference3 * VSL_low,
      pred_700_deaths_difference_VSL3_high = pred_700_deaths_difference3 * VSL_high,
      
      pred_550_deaths_difference_VSL1 = pred_550_deaths_difference1 * VSL,
      pred_550_deaths_difference_VSL1_low = pred_550_deaths_difference1 * VSL_low,
      pred_550_deaths_difference_VSL1_high = pred_550_deaths_difference1 * VSL_high,
      
      pred_550_deaths_difference_VSL2 = pred_550_deaths_difference2 * VSL,
      pred_550_deaths_difference_VSL2_low = pred_550_deaths_difference2 * VSL_low,
      pred_550_deaths_difference_VSL2_high = pred_550_deaths_difference2 * VSL_high,
      
      pred_550_deaths_difference_VSL3 = pred_550_deaths_difference3 * VSL,
      pred_550_deaths_difference_VSL3_low = pred_550_deaths_difference3 * VSL_low,
      pred_550_deaths_difference_VSL3_high = pred_550_deaths_difference3 * VSL_high
    ) %>%
    mutate(
      pred_900_deaths_difference_VSL1_discounted = pred_900_deaths_difference_VSL1 / (1 + 0.09) ^
        (Year - 2023),
      pred_900_deaths_difference_VSL1_low_discounted =  pred_900_deaths_difference_VSL1_low / (1 + 0.09) ^
        (Year - 2023),
      pred_900_deaths_difference_VSL1_high_discounted =  pred_900_deaths_difference_VSL1_high / (1 + 0.09) ^
        (Year - 2023),
      pred_900_deaths_difference_VSL2_discounted = pred_900_deaths_difference_VSL2 / (1 + 0.09) ^
        (Year - 2023),
      pred_900_deaths_difference_VSL2_low_discounted =  pred_900_deaths_difference_VSL2_low / (1 + 0.09) ^
        (Year - 2023),
      pred_900_deaths_difference_VSL2_high_discounted =  pred_900_deaths_difference_VSL2_high / (1 + 0.09) ^
        (Year - 2023),
      pred_900_deaths_difference_VSL3_discounted = pred_900_deaths_difference_VSL3 / (1 + 0.09) ^
        (Year - 2023),
      pred_900_deaths_difference_VSL3_low_discounted =  pred_900_deaths_difference_VSL3_low / (1 + 0.09) ^
        (Year - 2023),
      pred_900_deaths_difference_VSL3_high_discounted =  pred_900_deaths_difference_VSL3_high / (1 + 0.09) ^
        (Year - 2023),
      
      
      pred_700_deaths_difference_VSL1_discounted = pred_700_deaths_difference_VSL1 / (1 + 0.09) ^
        (Year - 2023),
      pred_700_deaths_difference_VSL1_low_discounted =  pred_700_deaths_difference_VSL1_low / (1 + 0.09) ^
        (Year - 2023),
      pred_700_deaths_difference_VSL1_high_discounted =  pred_700_deaths_difference_VSL1_high / (1 + 0.09) ^
        (Year - 2023),
      pred_700_deaths_difference_VSL2_discounted = pred_700_deaths_difference_VSL2 / (1 + 0.09) ^
        (Year - 2023),
      pred_700_deaths_difference_VSL2_low_discounted =  pred_700_deaths_difference_VSL2_low / (1 + 0.09) ^
        (Year - 2023),
      pred_700_deaths_difference_VSL2_high_discounted =  pred_700_deaths_difference_VSL2_high / (1 + 0.09) ^
        (Year - 2023),
      pred_700_deaths_difference_VSL3_discounted = pred_700_deaths_difference_VSL3 / (1 + 0.09) ^
        (Year - 2023),
      pred_700_deaths_difference_VSL3_low_discounted =  pred_700_deaths_difference_VSL3_low / (1 + 0.09) ^
        (Year - 2023),
      pred_700_deaths_difference_VSL3_high_discounted =  pred_700_deaths_difference_VSL3_high / (1 + 0.09) ^
        (Year - 2023),
      
      
      pred_550_deaths_difference_VSL1_discounted = pred_550_deaths_difference_VSL1 / (1 + 0.09) ^
        (Year - 2023),
      pred_550_deaths_difference_VSL1_low_discounted =  pred_550_deaths_difference_VSL1_low / (1 + 0.09) ^
        (Year - 2023),
      pred_550_deaths_difference_VSL1_high_discounted =  pred_550_deaths_difference_VSL1_high / (1 + 0.09) ^
        (Year - 2023),
      pred_550_deaths_difference_VSL2_discounted = pred_550_deaths_difference_VSL2 / (1 + 0.09) ^
        (Year - 2023),
      pred_550_deaths_difference_VSL2_low_discounted =  pred_550_deaths_difference_VSL2_low / (1 + 0.09) ^
        (Year - 2023),
      pred_550_deaths_difference_VSL2_high_discounted =  pred_550_deaths_difference_VSL2_high / (1 + 0.09) ^
        (Year - 2023),
      pred_550_deaths_difference_VSL3_discounted = pred_550_deaths_difference_VSL3 / (1 + 0.09) ^
        (Year - 2023),
      pred_550_deaths_difference_VSL3_low_discounted =  pred_550_deaths_difference_VSL3_low / (1 + 0.09) ^
        (Year - 2023),
      pred_550_deaths_difference_VSL3_high_discounted =  pred_550_deaths_difference_VSL3_high / (1 + 0.09) ^
        (Year - 2023)
      
    ) %>% 
    mutate(model = i)
  
  

  india_kg_pm_scenario_1_yr_long <- 
    india_kg_pm_scenario_1_yr_long %>% 
    bind_rows(india_kg_pm_scenario_1_yr)

  
  print(i)
  
  
}


write_rds(india_kg_pm_scenario_1_yr_long, "~/Downloads/india_kg_pm_scenarios_long.rds")





# BOOT PM2.5 KG response -------

for (i in 1:1000) {
  
  min1 = rnorm(1, 25, 5)
  max1 = rnorm(1, 100, 5)
  decline1 = rnorm(1, 70, 2.5)
  exp1 = rnorm(1, 4, .05)
  exp2 = rnorm(1, 8, .1)
  
  x=1:200
  PM25_1=min1 + (max1-min1) / (1 + (x/decline1)^exp1)
  PM25_2=min1 + (max1-min1) / (1 + (x/max1)^exp1)
  PM25_3=max1 - ((max1 - min1) / 200)*x
  
  kg_lpg_pm25 <- tibble(
    kg_lpg = 1:200,
    PM25_1 = PM25_1,
    PM25_2 = PM25_2,
    PM25_3 = PM25_3,
    model = i
  )
  
  kg_lpg_pm25_long <- 
    kg_lpg_pm25_long %>% 
    bind_rows(kg_lpg_pm25)

  print(i)
}


write_rds(kg_lpg_pm25_long, "~/Downloads/india_kg_lpg_pm25_long.rds")


